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Abstract 

In this paper we develop an algorithm to calculate the prices and Greeks 
of barrier options in a hyper-exponential additive model with piecewise con- 
stant parameters. We obtain an explicit semi-analytical expression for the first- 
passage probability. The solution rests on a randomization and an explicit matrix 
Wiener-Hopf factorization. Employing this result we derive explicit expressions 
for the Laplace-Fourier transforms of the prices and Greeks of barrier options. 
As a numerical illustration, the prices and Greeks of down-and-in digital and 
down-and-in call options are calculated for a set of parameters obtained by a si- 
multaneous calibration to Stoxx50E call options across strikes and four different 
maturities. By comparing the results with Monte-Carlo simulations, we show 
that the method is fast, accurate, and stable. 
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1 Introduction 



Barrier options are contracts whose pay-offs are activated or de-activated when the 
underlying process crosses a pre-specified level. These contracts are among the most 
popular path-dependent options. To value barrier options, a model needs to be 
sufficiently flexible to calibrate call option prices at different strikes and maturities. 
However, it is desirable to maintain a degree of analytical tractability to facilitate 
the calculations, especially for the Greeks or the sensitivities. These sensitivities 
describe the change in the model price with respect to a change in the underlying 
parameter, and are important for an appreciation of the robustness of the model's 
results. It is well known that the accurate evaluation of the Greeks is a challenging 
numerical problem, since standard PDE or Monte-Carlo methods are generally slow 
and unstable. 

It is well established that the geometric Brownian motion model lacks the flexi- 
bility to capture features in financial asset return data such as the skewness and the 
excess kurtosis. It cannot calibrate simultaneously to a set of call option prices. To 
address these limitations, one of the approaches consists of introducing jumps in the 
price process by replacing the Brownian motion by a Levy process. Levy models, 
such as the VG, CGMY, NIG, KoBoL, generalised hyperbolic, and Kou's double ex- 
ponential model, have been successfully applied to the valuation of European-type 
options. We refer to Cont and Tankov |14] . Boyarchenko and Levendorskii [7], and 
Schoutens [30] for background and references on the application of Levy models in 
option pricing. 

As observed by many authors, such as Eberlein and Kluge [16], or Carr and Wu 
, Levy models are generally not capable of calibrating option prices simultaneously 
across strikes and maturities. Empirical studies of S&P500 index data by Carr and 
Wu [11], and Pan [26], show that the implied jump intensities and the implied jump 
size distributions vary greatly over time. The prices of short-dated options exhibit 
a significantly larger risk-premium than that of long-dated options. This is reflected 
in the thicker tails of the implied marginal risk- neutral distributions, especially at 
short maturities. For example, in the equity markets, short-dated out-of-the money 
put options are relatively expensive since the risk of a large negative jump in the 
share is priced. Because of the stationarity and independence of the increments of a 
Levy process, the moments exhibit a rigid term structure that is different from what is 
observed in market data. This lack of flexibility can be overcome by considering mod- 
els driven by additive processes, which have independent and time-inhomogeneous 
increments. 

Additive models have been used for equity option pricing by Carr et al. [10], 
Galloway and Nolder j!7] , and by Eberlein and Kluge |16j for interest rate option 
pricing. Motivated by modelling considerations, Carr [10] proposed a self-similar 
additive model for the log-price, and reported good calibration results across time. 
Galloway and Nolder [IT] carried out a calibration study for various related models. 
Eberlein and Kluge [16] constructed an HJM model driven by an additive process 
with continuous characteristics, and they obtained a good fit for swaptions by using 
piecewise constant parameters. 
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In this paper we follow a similar approach: we model the share price by an ad- 
ditive process with hyper-exponential jumps. Hyper-exponential distributions are 
finite mean-mixtures of exponential distributions which can approximate monotone 
distribution arbitrarily closely. As first observed by Asmussen et al. [I], most of the 
popular Levy models used in mathematical finance possess completely monotone Levy 
densities and can therefore be approximated well by hyper-exponential Levy models. 
A hyper-exponential additive model is sufficiently flexible to allow for an accurate 
calibration to European option prices across strikes and multiple maturities. In ad- 
dition, if the parameters are piecewise constant, the model admits semi-analytical 
expressions for prices and Greeks of barrier options. 

There currently is a body of literature devoted to various aspects of pricing barrier 
options. In the setting of Levy models, a transform-based approach to price barrier 
options has been developed in a number of papers, including Geman and Yor |18j . 
Kou and Wang [21 1, Davydov and Linetsky [15] . Boyarchenko and Levendorskii [8]. 
In particular, Kou and Wang [21], Kou et al. [22], Sepp [29], Lipton [21], and Jeannin 
and Pistorius pj5] considered the cases of Levy processes with double-exponential and 
hyper-exponential jumps. 

In this paper, the transform algorithm that we develop is based on a so-called 
matrix Wiener-Hopf factorization. Such matrix factorizations were first studied by 
London et al. [25] and Rogers [27] for (noisy) fluid models. Jiang and Pistorius 
[20] developed matrix- Wiener factorization results for regime-switching models with 
jumps. We show that by suitably randomizing the parameters the distributions of the 
infimum and supremum of the randomized hyper-exponential additive process can be 
explicitly expressed in terms of a matrix Wiener-Hopf factorization. We use these 
results to derive semi-analytical expressions for the first-passage time probabilities, for 
the prices, and for the Greeks of barrier options, up to a multi-dimensional transform. 
The actual prices are subsequently obtained by inverting this transform. 

As a numerical illustration, we calibrate the hyper-exponential additive model 
to Eurostoxx prices quoted on 27 February 2007 at four different maturities. We 
calculate in this setting down-and-in digital and down-and-in call option prices and 
Greeks (delta and gamma). To invert the transform, we use a contour deformation 
algorithm and a fractional Fast Fourier Transform algorithm, developed by Talbot 
[5T] , Bailey and Swarztrauber [6], and Chourdakis [12], [13] ■ We also compare it to 
Monte-Carlo Euler scheme simulations. We find that the algorithm is accurate and 
stable, and much faster than Monte-Carlo simulations (especially for the Greeks). 
This method is suitable for applications in which the number of periods is not too large 
(up to four). When a larger number of periods is required, the direct inversion method 
used here is no longer feasible. The subject still needs to be further investigated and 
is left for future research. 

The remainder of the paper is organized as follows. In Section [2] we define the 
hyper-exponential additive model and present its application to European call option 
pricing. In Sections [3] and H] we derive semi-analytical expressions for the first-passage 
probabilities of a hyper-exponential additive process in terms of a matrix Wiener-Hopf 
factorisation, and for the prices and Greeks of barrier options. In Section[S]we present 
numerical results. 
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2 The model 



2.1 Additive processes 

We consider an asset price process S modelled as the exponential 

St = S e Xt 

of an additive process X. Informally, an additive process can be described as a Levy 
process with time-dependent characteristics or, equivalently, as a process with inde- 
pendent but non-stationary increments. We briefly review below some key properties 
of additive processes. For further background on additive processes and their appli- 
cations in finance, we refer to Sato [28], and to Cont and Tankov |14j . An additive 
process can be defined more formally as follows. 

Definition 1 For a given T > 0, X = {X t ,t e [0, T]} is an additive process if 

(i) X = 0, 

(ii) For any finite partition < to < t\ ■ ■ ■ < t^ < T , the random variables Xt h — 
X tk _ 1 , ■ ■ • ,X ta are independent, 

(Hi) The sample paths t t-> Xt have cddldg modifications almost surely. 

If X is an additive process, then, for every t £ [0, T], Xt has an infinitely divisible 
distribution with Levy triplet (Mt, T^,At); that is, the characteristic function of Xt 
is given by <&t(u) = exp^^-u) ]. According to the Levy-Khintchine formula, is the 
characteristic exponent given by 

* t (u) = iuM t - -j-v 2 + J [e lux - 1 - mxl {W<1} j A t (dx), 

with Mt, Tit £ and where At the Levy measure satisfies the integrability constraint 

J (1 A x 2 )A t (dx) < oo. 

The law of the additive process {X t ,t £ [0, T]} is determined by the collection of 
Levy triplets {(M t ,S|,A t ) for t £ [0,T]}. If the Levy triplets are time- independent, 
A is a Levy process. If the additive process has absolutely continuous characteristics, 
the Levy triplets take the explicit form 

M(t) = [ n(s)ds, E(t) = f a 2 (s)ds, 
Jo Jo 

At(B) = / / g(s,x)dxds for Borel sets B, 
Jo Jb 

where /x, a 2 : [0, T] — > M and g : [0,T] xl->l are integrable functions, with g and 
cr 2 non-negative. We call the functions (fi,a 2 ,g) the local triplet of X. 
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We assume that we have been given deterministic integrable functions r(t) and d(t) 
representing the short rate and the dividend yield, and that the characteristic expo- 
nent of Xt satisfies 

*t(-i)= f [r(s)-d(s)]ds, (1) 
J o 

or equivalently 

a 2 (t) f°° 

+ "4^ + - 1 - l W< iyx]g(t, x)dx = r{t) - d{t). (2) 



2 

^ J — oo 

It follows that the discounted process 

5" exp( / \r{s) - d(s)]ds) 
Jo 

is a martingale if and only if ([!]) (or, equivalently, ([2])) holds. 



2.2 Hyper-exponential additive processes 

In what follows we restrict the discussion to a hyper-exponential additive process X 
which is specified by its local triplet (fi, a 2 ,g) where g is given by 

n + n~ 

g(t,x) = ^vr+(t)a+(t)e-^Wn {l . >0} + ^vr7(t)a-(t)e-^WI-ll {x<0} , 

k=l j=l 

where vr^(t) and arr(t) are non- negative. The continuous part of X consists of a 
diffusion with time-dependent drift fx(t) and volatility cr(t). The jump part of X is of 
finite activity and forms an inhomogeneous compound Poisson process where positive 
and negative jumps occur at the rates 

n+ n~ 

\ + {t) := 5>+(*) and A"(t) := J^^C*). 

k=l j=l 

and jump sizes are distributed according to a hyper-exponential distribution. 
Small random price movements are intuitively modelled by the diffusion part, whereas 
sudden changes of the price are captured by the jump-part of X. If we take n ± = 1, 
the jump-sizes are exponentially distributed, and this model reduces to an extension 
of the Kou model with time-dependent parameters. 



2.3 Piecewise constant parameters 

To reduce the dimension of the available parameter set, we take the functions fi(t), cr(t) 
and g(t, •) to be piecewise constant. Given that we have a finite set of European call 
options with different maturities Ti, . . . , 2jv 3 we take the local parameters to be con- 
stant between the different maturities Tj. Then for all t £ (Tj_i,Tj], (with To = 0) 
we set 

Kt) = ^\ a 2 (t) = a 2 ^\ g(t,x) = g^(x), i = l,...,N. (3) 
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For t € (Tj_i,Tj] the characteristic exponent of X t — Xx i _ 1 is given by 

Vt^M =■ * (i) («)> 

where 

*<•>(„> = ,<•)„. - ^ + y>:«> (^L-) ( . (4 ) 

^ k=l \ a k ~ U1 J 3=1 \ a i +W 

3 First passage probabilities 

The value of a digital barrier option can be expressed in terms of the distribution 

F(+)(x;T) = P(X T < x) 

of the running supremum 

Xt = sup X s 

s<T 

of X, or equivalently, the distribution of the first-passage time 

T + (x) = inf{t > : X t > x} 

which is related to by 

P(T + (x) < T) = 1 - F {+ \x; T). 

Whereas for a Levy process the distributions of the infimum and supremum are linked 
to the characteristic exponent by the so-called Wiener-Hopf factorization, such a 
result does not exist for general additive processes, because of the time-dependence 
of the parameters. However, in the case of piecewise constant parameters, the triplet 
changes only at deterministic times, so that as a consequence the distribution function 
F(+)(x) = TW, . . . , fW) only depends on the inter-jump times TW = T; - 

T-i (with To = 0). In this case, as we show below, the iV-dimensional Laplace 
transform G^ + \x,q) of , given by 

where q = (qi, . . . ,qN), is expressed explicitly in terms of a matrix Wiener-Hopf 
factorization. To state this result we need to introduce some further notation. 

For any vector v = (v±, . . . ,v n ), we denote by A„ the diagonal matrix A v = 
(vi, i = l,..., n)diag- Let Q be the N(l + n + + n~) x N(l + n + + n~) matrix given 
in block notation by 



where 

t+ T+ 



^ + = ( g ; Aa £V (6) 
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Here A = (A+ + A: 



1, . . . , N), and G and b + are the N x N and N x Nn^ 



matrices in block notation given by 



G 



f-Ql Qi 

-Q2 Q2 



\ 



—Qn-i Qn-i 
-Qn/ 



/yKU 



b + 



7T 



+ (2) 



V 



7T 



• (7) 



Here 7r + W is the row- vector 7r + W = (7r, + , 1 = 1,. 
by 



, n 



I, and where t + , T + are given 



/-A, 



V 



\ 



v 



(8) 



where cr is the column-vector a 4 

c- = (r o-), 



(af,i = 1, . . . , n + )', and 
6" 



+ ; ' 



(9) 



where are n ± A r x n ± A r zero matrices, and b~ , T~ , and t~ are given by (J7J) and 
([8j) with 7T + ^^ and a + replaced by 7r~^ and a~ . The matrix Q is a generator 
matrix, that is, a square matrix with non-negative off-diagonal elements and non- 
positive row sums, and defines a Markov chain. This Markov chain is associated to 
a randomization and embedding of the additive process X (which will be illustrated 
with a concrete example below). We recall that a sub-probability matrix is a matrix 
with non-negative elements and row sums not larger than one. By applying the 
matrix Wiener-Hopf factorization results of Jiang and Pistorius |2U] to the current 
setting we arrive at the following conclusion. 



Theorem 1 It holds that 
where 

ei = (1,0,..., 0) 



1 



qi . . . q N 



[1 



/ a Q+x 



e-i e 



1 , 



(10) 



and 



! = (!,.. .,!)'. 



is an N(l + n + ) x iV(l + n + ) generator matrix that together with r] + an Nn x 
N(l + n + ) sub-probability matrix, solves the system of matrix equations 



( \S 2 Q\ - V+Q+ + H+ + D-rf* 



O, 



(11) 



-r, + Q + + C- +T-rf 



O. 
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Here the O 's are zero matrices of appropriate sizes, and in block notation we have, 



S 2 



A, 



O 



+ ; ' 



V 



r+ ' 



(12) 



with 



a- = {a 2 ,i = 1,...,N), p = (n h i = 1,...,JV). 

+ and I + represent n + N x n + N zero and identity matrices, respectively. 

By applying Theorem [T] to —X, we find the corresponding pair of matrices (Q_,?7_). 
The quadruple (Q + ,rj + , Q-,r)J) is called a matrix Wiener-Hopf factorization of Q. 

Example. To illustrate this approach, we consider a hyper-exponential additive 
process X on [0, T2] whose parameters are constant during the periods [0, T\] and 
[Ti,T2]. In the first period X evolves as a jump- diffusion with positive and nega- 
tive exponential jumps with means and jump rates l/a + , A + and l/a~,A~. In the 
second period X is a Brownian motion with drift. The idea is to randomize the 
times between maturities by replacing = T\ and = T2 — T\ with inde- 
pendent exponential random variables having means q^ 1 and q^ 1 - This results in a 
regime-switching jump-diffusion with the regime only jumping from state 1 to state 
2, according to the generator matrix 



G 



-<?i 




-Q2 



We associate to the regime-switching process a continuous Markov additive process, 
which can be informally obtained by replacing positive and negative jumps with 
stretched slopes of +1 and —1 (see Asmussen [3] for background on this embedding). 
As described in [20], in this case the generator of the modulating Markov process is 
given by 



Q 



H+ D- 
C~ T~ 



( 


- qi -x+-x- 


qi 


A+ 


A~ 









-Q2 












a+ 





-a+ 







\ 










— a~ 


) 



with the matrices S 2 and V + in Theorem [T] given by 



07 



S 2 



Mi 



V 



f'2 



3.1 Solution of the matrix equation 

To solve the system (jlip . which is a Ricatti-type matrix equation, we follow a spectral 
approach and determine the spectral decomposition of Q + . Denoting by h{p) a 
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Figure 1: Paths of the various processes related to log-price process X are illustrated. 
Here X is a hyper-exponential additive process on the period [0, T2] whose param- 
eters are constant during the periods [0,Ti] and [Ti,T2]. During the first period 
[0,Ti], the process X evolves as a jump- diffusion with volatility a\, drift fi± = and 
exponentially-distributed jumps. The positive and negative jumps are exponentially- 
distributed with means and jump rates l/a + , A + and 1/oT , X", respectively. During 
the second period [T\ , T2] , X evolves as a linear Brownian motion with volatility 02 
and drift \xi = 0. Associated to X is a continuous Markov additive process A, which 
can be obtained from X by replacing the positive and negative jumps with linear 
stretches of path of slopes +1 and —1, and replacing the fixed times T%, T2 — T\ by 
independent exponential random times ei and &2 with parameters qi,q2- The process 
Y that records the current state or regime of A is a Markov process with generating 
matrix Q. When Y takes values 1 and 2, A evolves as a linear Brownian motion 
with zero drift and volatility o\ and 02 respectively; and when Y is 3 and 4, A is a 
positive or negative unit drift. These linear stretches of paths of A originate from 
the jumps of X. A jump of Y from one state to another is induced either by a jump 
of X or by a switch of the set of parameters that determine the dynamics of X. By 
time-changing A by the time To(t) up to time t that Y was equal to 1 and 2, we 
recover a regime-switching jump-diffusion; that is, the process {A(To(t)),t > 0} is in 
law equal to a regime-switching jump-diffusion. Finally, replacement of the times at 
which a regime switch occurs by T± and T2 results in a process that has the same 
distribution as X. 
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(column) eigenvector of Q + corresponding to the eigenvalue p, one finds that it is a 
matter of algebra to verify that the system (jlip can be equivalently rewritten as 

H; + )q 2 + -K; + H +q C')= o - 

Here O is a N(l + n + + n~) square zero matrix, I is an N(l + n + ) identity matrix, 
and 

-v 2 - I o- I . v - f /- 



Defining the matrix iT(s) by 



s 2 

K(s) = —S 2 + sV + Q, (13) 
we find that h(p) solves the linear system 

K[- P ] f+) Hp) = o, 

which implies that p is a root of the equation det K(s) = 0. The following result 
characterizes the eigenvalues of Q + (see Appendix A): 

Lemma 1 

(i) It holds that 

N 



|det(tf( fl ))|=n 



1=1 



|*W(-b) - «| J] |s - a+| J] \s + af| \, (14) 



k=l 1=1 



where jg given in (J3|). 
fiij T7ie equation 

detiq» = (15) 
/ias JV(1 + positive roots and N(l + re - ) negative roots. 

Since — Q + is the negative of a generator matrix, it is non-negative definite, so its 
eigenvalues are non-negative and are given by the positive roots of (|15|) . In particular, 
if the positive roots 

p + = {p+:i = l,...,N(n+ + l)) 
of equation (|15j) are distinct, it follows from Lemma [1] and Theorem [1] that 

G (+) (x,q) = (qi-'-qN)- 1 x [1 - ejtf+e"^^ 1 !], 
where C/+ = (h(pf),i = l,...,N(n + + 1)). 
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3.2 The final position and the first exit time 

The valuation of barrier options involves the joint distribution of the final position at 
maturity T and the first exit time. We will extend the results in the previous section 
by considering the following: 

= E[e sXT l {T+{x)<T} ]. 

F^ + \x,s) depends on time only through the inter-maturity times (T^\ . . . ,T^). 
The Laplace transform G^ + \x, s; q) of F^ + \x, s) in (T^\ . . . , T^), can be expressed 
in terms of Q + and K{s) as follows: 

Proposition 1 It holds that 

G {+) (x, s, q) = x e' 1 e Q+x i^(s)- 1 ^(0)l (16) 

Ql • • • QN 

for all s € C with Re(s) G (— min J=1 „- a~,min fc=lj n + ac£). 
A proof is given in Appendix A. 

3.3 First passage to a lower level 

The form of the analogous distributions concerning the infimum 

F(-\x) = P(-X T <x), F { -\x,s) = E[e sXT l { _x T>x} ], x > 0, 

can be found by applying the results in the previous section to the process —X. More 
specifically, it is straightforward to check that the A^-dimensional Laplace transforms 
G(~\x, q) and \x, s, q) are given by (fTU|) and (|16|) replacing Q + by Q~ . (Q~,rj~) 
satisfies the system of matrix equations (jlip with V + , H + , D~ ,C~ ,T~ replaced by 
V~ , H~, D + , C + , T + , where the latter set is defined by interchanging + and — in 
equations (jI2]> . ([9]), dHJ) and (j6j). It is straightforward to verify that (i) an eigenvector 
h(p) of Q~ corresponding to eigenvalue p satisfies 

K[ P ] Hp) = o, 

where I is an N(l + n~) identity matrix, and (ii) that, in view of Lemma [TJ the 
eigenvalues of Q~ are given by the negative roots 

P = (pJ,j = l,...,N(pr + lf) 

of det K{p) = 0. 
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4 Prices and Greeks of digital and barrier options 



Using the first-passage results from the previous section we derive semi-analytical 
expressions for the prices and sensitivities of down-and-in digital and knock-in call 
options. A down-and-in digital option at level H < So is a contract that pays out one 
unit at maturity T if the price S has down-crossed the level H before T. Similarly, a 
down-and-in call option at level H < So and with strike K is a call option whose pay- 
off is activated once S down-crosses H. Taking the risk-free rate r and the dividend 
rate d to be constant, the arbitrage free prices of a down-and-in digital (DID) and 
a call option (DIC) are given respectively by 

DID(T,H,S ) = e~^ T E [l {inf ^ T Ss<H} ] = e~^ T P(X T < h), 

where h = \og(H / S$) is the log-barrier, and 

DIC(T,H,K,S ) = e~^ T S E [(e XT - e k )+l { x T<h} ] , 

where k = log(K/So) denotes the log-strike. Let DID(q) denote the joint Laplace 
transform of DID in the inter-maturity times (T^\ . . . , T^ N ') (with T( N > = T), and 
denote by DIC (q, s) the Laplace- Fourier transform in (T^\ . . . ,T^ N ^) and in the 
log-strike k. Then we have the following result: 

Proposition 2 For h = log (I?/ So) < it holds that 

DID(q) = -j—xe\e Q ~ h l, (17) 

DIC(q,s) = ^ x ^ e q" fc Jf (0-^(0)1, (18) 

where k = \og(K/So), q = (gi, . . . ,qx), andb = a+is+1, c(q) = (qx+r) ■ ■ ■ (qN + r). 

Before we give the proof we observe that from the explicit expressions (|17p and 
(I18D semi-analytical formulas can be obtained for the delta and gamma of the down- 
and-in digital and call options (i.e. the first and second derivatives of the option value 
with respect to the spot So). Indeed, the derivatives of the expressions (fTT|) and (|18p 
with respect to So are equal to the Laplace-Fourier transforms of the derivatives of 
the option, as integration and differentiation are interchangeable in this case. In the 
case of a down-and-in digital option we find that the Laplace transforms £±did and 
Tdid of the delta &did and gamma Tdid are given by 

A DID (d) = — -L- x e[Q-e^ h l, f DID (q) = — x e[[(Q-) 2 + Q~]e Q ' h l. 
c{q)S c{q)S^ 

Proof of Proposition \B[ The expression (|17|) is a direct consequence of Theorem [T] 
(see also Section [3.3p . To verify (|18p we start by taking the Fourier transform in k 
and find as in (|2ip that the Fourier transform DIC* is given by 

D/C*(a + i S ) = F ^ ( ~y ) , (19) 
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where b = a + is + 1 and 



F y >(x,b) = E[e bX ^l { _x T > x }}- 



From Proposition [T] we deduce that the form the joint Laplace transform G(x, b, q) 



5 Numerical results 
5.1 Calibration 

To determine a parameter set to test the method, we calibrate the hyper-exponential 
additive model to Eurostoxx call options at four different maturities, observed in 
the market on 20 February 2007. The spot price is EUR 4150, the risk-free rate is 
assumed to be fixed at r = 0.03, and the dividend rate is taken to be zero. As we 
find that inclusion of positive jumps does not substantially improve the calibration 
results, we only consider negative jumps, and we specify the jump size parameters 
to be (oti , (X2) = (3, 10). The jump arrival rates tv~ and the volatility are piecewise 
constant in time, and are estimated by minimizing the root-mean-square error be- 
tween model and observed market call prices. Using the well known Fourier transform 
method (briefly recalled in Appendix [B]) the calibration is carried out maturity by 
maturity under constraints through a bootstrapping method with well-defined local 
triplets: 

i Calibrate call prices at T\ to obtain the parameters (a(Ti), it^(Ti), n^(T\)). 

ii For j = 2, . . . , N calibrate call prices at Tj to obtain (a(Tj),TT^(Tj), ir^(Tj)). 

In Figure 2 the calibration results are presented with plots of the market and 
model implied volatility surfaces corresponding to the four maturities 6m, 1Y, 3Y 
and 5Y. The root-mean-square error (RMSE) and the average relative percentage 
error (ARPE) are equal to 5.30 and 1.1%. We compare it to a price process that 
follows a Levy process with hyper-exponential jumps (i.e. with constant parameters 
over time), and find that the calibration of the four maturities in that case give a 
RMSE of 9.82 and an ARPE of 2.9%. In Table 1 the resulting parameter sets are 
displayed under the hyper-exponential additive and Levy models. In the case of the 
hyper-exponential model we observe a high jump intensity of small jumps for short 
maturities that decrease substantially over time. This is consistent with the finding 
of Carr and Wu [llj, and Pan [26]. 




x e' 1 e c? " /l K(s)- 1 J ft:(0)l. 



(20) 



Combining (fT9|) and (f20|) completes the proof. 



□ 
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Figure 2: Calibration of Eurostoxx call option prices on February 20th 2007 for 
an additive hyper-exponential model with piecewise constant parameters. Crosses 
represent market implied volatilities and circles model implied volatilities for four 
maturities: 6 months (red), 1 year (orange), 3 years (green), 5 years (blue). 
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Additive hyperexponential model 




Tstart 


TEnd 


a 




^2 





0.5 


0.0995 


0.0371 


11.1819 


0.5 


1 


0.0759 


0.2091 


9.9540 


1 


3 


0.0786 


0.4738 


7.0322 


3 


5 


0.0858 


0.8084 


0.2361 


Levy hyperexponential model 




Tstart 


TEnd 


a 




^2 





5 


0.1171 


0.5693 


0.0165 



Table 1: Calibration of Eurostoxx call option prices quoted on February 20 2007 for 
an additive hyper-exponential model with piecewise constant parameters and for a 
Levy hyper-exponential model with jump parameters = 3 and = 10. The 
interest and dividend rates are r = 0.03 and d = 0. The maturities are given in 
years. 
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5.2 Results for the barrier option prices and Greeks 

Using the parameter set found in the calibration of the Eurostoxx call options, we 
value barrier and digital options on the Eurostoxx index, modelling its price pro- 
cess as the exponential of a hyper-exponential additive process. We use the semi- 
analytical results in Proposition [2j To invert the multi-dimensional Laplace trans- 
forms we choose Talbot's method |31j (see also |12j ) for down-and-in digital options. 
We combine it with the fractional FFT algorithm of Bailey and Swarztrauber [BJ, 
and Chourdakis [T3] for down-and-in call options. See Appendix O for a detailed 
description of the implementation of these transform algorithms. We compare it to 
the same quantities calculated by Monte-Carlo simulations, using a standard Euler 
scheme. 

5.2.1 Down-and-in digital options 

We price down-and-in digital options with a maturity of five years for different spot 
levels. We evaluate the required 4-dimensional Laplace transform over two time 
increments of six months, and two of two years, that is, = 0.5, = 0.5, 

T (3) _ 2 5 r( 4 ) = 2. We use Talbot's algorithm with M = 6 (see Appendix O for 
an explanation of this parameter). Using Mathematica to run the algorithm, the 
computation time was five minutes on a 3189 Mhz computer to calculate prices and 
Greeks for fourteen different spot levels. The calculation of first passage probabilities 
using Monte-Carlo simulations requires a large number of time steps and paths. We 
use one million paths with St = 5 x 10 -5 and it takes several hours to obtain stable 
Greeks in C++. Error bounds cannot be obtained analytically, but we observe in 
Table [2] that the results of the transform method agree with Monte-Carlo simulation 
results. Figures 3 and 4 report prices and Greeks for down-and-in digital options. The 
options are expressed as a percentage of the spot price. The values of the sensitivities 
are expressed as fractions of the spot price So- 

This transform algorithm is particulary efficient at a book level, since once the 
generating matrices Q~ of the infimum have been calculated for different values of 
the vector q the calculation of prices and Greeks of any digital barrier product is just 
a matter of summation. 

5.2.2 Down-and-in call options 

We value down-and-in call options with a maturity of one year for different strike 
levels. In this two-dimensional Laplace inversion is required over time incre- 

ments TW = 0.5 and = 0.5. For the inversions of the Laplace transform and 
the Fourier transform, we set M = 7 and N = 1024 (refer to Appendix O for an 
explanation of these parameters). For Monte-Carlo simulations, we use one million 
paths with time step St = 2.5 x 10 -5 . The option prices and Greeks obtained by 
the two methods are reported in table [3] and figures 0] and We observe that the 
results of the transform method agree with the Monte-Carlo simulation results. Us- 
ing Mathematica again to run the algorithm, the computation time is ten minutes to 
calculate the option prices, delta and gamma for eleven different levels of the strike. 
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Since the option prices and Greeks of a down-and-in call option are obtained via a 
Fourier-Laplace transform, it takes more time than in the case of a digital option 
(approximately twice as long), which is still much faster than a Monte-Carlo Euler 
scheme. We note that the transform algorithm is particulary efficient for the pricing 
of options with different strikes, as we obtain by FrFFT inversion the prices and 
Greeks of any down-and-in call options on a log-strike grid. 
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Figure 3: Prices of down-and-in digital options with a barrier H set at 90% of EUR 4150 
and maturity T — 5 years. Semi-analytical results are indicated with the symbol x and 
Monte-Carlo results with the symbol o. 
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Figure 4: Greeks of down-and-in digital options with a barrier H set at 90% of EUR 4150 
and maturity T = 5 years. Semi-analytical results are indicated with the symbol x and 
Monte-Carlo results with the symbol o. 
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Figure 5: Prices of down-and-in call options with a barrier H set at 90%, a spot at EUR 
4150 and maturity T = 1 year. Semi-analytical results are indicated with the symbol x and 
Monte-Carlo results with the symbol o. 
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Figure 6: Greeks of down-and-in call options with a barrier H set at 90%, a spot at EUR 
4150 and maturity T = 1 year. Semi-analytical results are indicated with the symbol x and 
Monte-Carlo results with the symbol o. 



17 





Down-and-in digital options 




Price 


Delta (xlO~ J ) 


Gamma 


(Xl0~ 7 ) 


% 




TA 


(MC- 


, MC+) 


TA 


MC 


TA 


MC 


92 





7261 


( 0.7000, 


0.7499) 


-1.226 


-1.232 


7.34 


7.54 


94 





6446 


( 0.6137, 


0.6735) 


-0.812 


-0.819 


3.64 


3.70 


96 





5893 


( 0.5566, 


0.6207) 


-0.577 


-0.579 


1.92 


1.94 


98 





5478 


( 0.5143, 


0.5805) 


-0.450 


-0.451 


1.09 


1.11 


100 





5140 


( 0.4800, 


0.5475) 


-0.374 


-0.377 


0.67 


0.67 


102 





4850 


( 0.4506, 


0.5189) 


-0.328 


-0.330 


0.46 


0.45 


104 





4592 


( 0.4245, 


0.4932) 


-0.294 


-0.297 


0.34 


0.33 


106 





4358 


( 0.4008, 


0.4697) 


-0.269 


-0.271 


0.28 


0.30 


108 





4144 


( 0.3795, 


0.4483) 


-0.247 


-0.247 


0.24 


0.24 


110 





3946 


( 0.3599, 


0.4285) 


-0.229 


-0.228 


0.21 


0.21 


112 





3762 


( 0.3518, 


0.4010) 


-0.212 


-0.212 


0.19 


0.19 


114 





3592 


( 0.3250, 


0.3917) 


-0.198 


-0.197 


0.17 


0.17 


116 





3433 


( 0.3095, 


0.3771) 


-0.184 


-0.182 


0.15 


0.16 


118 





3285 


( 0.3012, 


0.3632) 


-0.172 


-0.170 


0.14 


0.14 



Table 2: Down-and-in digital options with barrier H set at 90% of EUR 4150. The first 
column contains the spot price as a percentage figure of 4150. The columns with TA contain 
the results obtained using the transform algorithm, whereas MC refers to Monte-Carlo results. 
In the case of the price the Monte-Carlo results are reported in the form of a 95% confidence 
interval, and in the other cases as a point estimate. 







Down-and-i 


n call opt 


ons 








Price 


Delta (X10 


- 1 ) 


Gamma 


(X10~ 4 ) 


% 


TA 


(MC- , MC+) 




TA 




MC 


TA 


MC 


80 


111.13 


( 108.85, 111.20) 


-2 


458 


-2 


440 


9.49 


9.39 


82 


93.48 


( 91.46, 93.55) 


-2 


124 


-2 


107 


8.43 


8.34 


84 


77.14 


( 75.35, 77.19) 


-1 


808 


-1 


792 


7.45 


7.37 


86 


62.28 


( 60.71, 62.33) 


-1 


514 


-1 


499 


6.47 


6.40 


88 


49.08 


( 47.71, 49.13) 


-1 


242 


-1 


231 


5.57 


5.51 


90 


37.63 


( 36.7, 37.80) 


-0 


991 


-0 


988 


4.65 


4.59 


92 


28.11 


( 27.43, 28.46) 


-0 


783 


-0 


776 


3.88 


3.84 


94 


20.94 


( 20.13, 21.01) 


-0 


599 


-0 


593 


3.07 


3.03 


96 


15.11 


( 14.45, 15.16) 


-0 


448 


-0 


443 


2.45 


2.43 


98 


10.66 


( 10.12, 10.72) 


-0 


326 


-0 


320 


1.85 


1.84 


100 


7.36 


( 6.92, 7.41) 


-0 


231 


-0 


226 


1.37 


1.36 


102 


4.97 


( 4.63, 5.02) 


-0 


159 


-0 


156 


0.94 


0.93 


104 


3.28 


( 3.03, 3.34) 


-0 


107 


-0 


104 


0.64 


0.63 


106 


2.12 


( 1.94, 2.18) 


-0 


070 


-0 


069 


0.43 


0.42 


108 


1.35 


( 1.22, 1.41) 


-0 


045 


-0 


044 


0.28 


0.28 


110 


0.84 


( 0.74, 0.89) 


-0 


028 


-0 


027 


0.23 


0.22 


112 


0.51 


( 0.47, 0.54) 


-0 


017 


-0 


017 


0.18 


0.18 


114 


0.30 


( 0.25, 0.34) 


-0 


010 


-0 


009 


0.13 


0.13 


116 


0.18 


( 0.16, 0.19) 


-0 


006 


-0 


006 


0.10 


0.10 


118 


0.10 


( 0.08, 0.11) 


-0 


003 


-0 


003 


0.07 


0.07 


120 


0.06 


( 0.03, 0.06) 


-0 


002 


-0 


002 


0.03 


0.03 



Table 3: Down-and-in call options with a barrier H set at 90%, a spot at EUR 4150 and 
maturity T = 1 year for a range of strikes. The first column contains the strike level as a 
percentage figure of the spot EUR 4150. The columns with TA contain the results obtained 
using the transform algorithm, whereas MC refers to Monte-Carlo results. In the case of the 
price the Monte-Carlo results are reported in the form of a 95% confidence interval, and in 
the other cases as a point estimate. 
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APPENDIX 



A Proofs 

Proof of Lemma d' 

(i) It is straightforward to verify that K#(s) can be obtained from K(s) by inter- 
changing some columns and rows, where K^(s) is given by 



K#(s 



( K x (s) D a 

K 2 (s) D 2 



\ 



\ 



K N (s) J 



where K w (s) and D w are square matrices of dimension n + + n +1. There are 
defined respectively by 



K w ( 8 ) 



a-. 



V 



a 



TT 1 (w) 




7T m {w) 


irl (w) 




7T+(U> 


—s — 








































S Qt-rn 




















s — af 








































s — a 



where A J = Yli^fiw), and 



(D 



q w if i = j = 1 
otherwise 



Therefore | det(-fC(s))| is equal to | det(K#(s)\. To proceed we recall an identity 
from matrix algebra. Let M be a matrix of the form 



M 



Mn M12 
M 2 i M 22 



in block notation, where M22 is invertible. Then 

det(M) = det(M 22 ) det(Mn + M X2 M^ M 2 i). 
Using this identity, it is a matter of algebra to verify by induction that 
det(iT # (s)) = det(Ki(s)) • • • det(K N (s)). 
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As a consequence we find, by applying this matrix identity, that 



2 n + n q/ — 

det(K w (s))= I ^ 2 + ^ + J>+M^^ + J>7W ~ " A - 



H" 1 



x Y[(s - a+) JJ(-s - ^ 

i=l j=i 



and the assertion follows in view of (HJ). 



(ii) Using the intermediate value theorem and the specific form of it is straight- 
forward to check that the equation ui) = q, q > has n + 1 positive roots 
and m + 1 negative roots , satisfying 

/Wl < ~ Q m < Pm < • • • < -«1 < Pi < < ^ + < a l < ■ • • < Pn < a n < Pn+V 

Since det(K w (s)) is a polynomial of degree n+m+2, it follows that all the roots 
of det(K w (s)) = are given by (pf, i = 1, . . . , n + 1) and (pj ,j = 1, . . . , m + 1). 
In view of the form of det(K(s)) derived in (i) the assertion follows. 

□ 



Proof of Proposition^ Consider the following randomization of X obtained by ran- 
domizing the inter-maturity times by replacing them by independent exponential 
random times with means q7 , and call this process X. The process X is a regime- 
switching jump-diffusion, where the only regime switches that can occur are from i 
to i + 1 at rate qi (i = 1, ... ,N — 1), and from the final state N to an absorbing 
'graveyard state' d. As shown in (20], the process X is equal to a time-changed con- 
tinuous process A, say. Denoting by ( the epoch at which A is sent to d, by V the 
modulating Markov chain, and by r = inf{£ > : A T = x}, we have 

qi ---q N G(x,s,q) = E[e sA <-l { -j ( _ >x} ] 
= E[e sA <-l {T<} ] 
= E[e sA ^l {T<0 f(Y T )] 
= e sx E[l {T<c} f(Y T )], 

with 

f(y) = E[e sA t-\A = 0,Y = y], 

where the last two lines follow by the Markov property of (A, Y) and the fact that A is 
continuous. To guarantee that all the expressions are well defined in this calculation 
s has to be such that i?[e *] < oo, which corresponds to the restriction that 

Re(s) G (— min a~ , mina+). 

3 i 

In |20j it was shown that the vector f = (f(y), y £ N), where N denotes the state 
space of Y, is given by 

f = K(s)- 1 Ql, 
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where the matrix K(s) is given in ()13p . Combining these results with Theorem [T] we 
find that 

qx ■ ■ ■ q N G(x, s, q) = e[e Q+x K^Ql, 
and the proof is complete. □ 

B European call options 

Under the hyper-exponential additive model with piecewise constant parameters Q , 
the characteristic function at time T is explicitly given by 



= exp 

with tfCi) as given in (j3J). The price of a European call with maturity Tj can thus 
be efficiently calculated using a well-established Fourier transform method, which we 
briefly recall. The Fourier transform C™. over k of C^^k), the price of a call option 
with log-strike k = log(K / Sq) and maturity Tj, can be explicitly expressed in terms 
of the characteristic function $W(u) as follows: 

/oo 
-oo 

= ■ ,*"'!» T'° t 1 '!' . (21) 

(a + iu)(a + 1 + iu) v 7 

Since the call pay-off function itself is not square-integrable in the log-strike, the axis 
of integration is here shifted over ia which corresponds to exponentially dampening 
the pay-off function at a rate a, which is usually taken to be a = 0.75 (see Carr 
and Madan [9]). The call option prices are then determined by inverting the Fourier 
transform: 

C Ti (k) = ^ / e~ lfc V w T } .\ dv. (22) 

lV 7 tt Jo (a + iv)(a + l + iv) V ; 

C Transform inversion algorithms 
C.l Multi-dimensional Laplace inversion 

To evaluate down-and-in digital option prices (DID), we invert the multi-dimensional 
Laplace transform ()17j) to obtain 

DJZ)(5 0j /» J T) = ^p jT ■■■ j c e^ + - + ^ TN DID(S ,h,ci)dq. (23) 

where T = (T\, . . . , Tjy) and C ra are vertical lines in the complex plane defined by 
q-n = f n + iy n for n = 1, . . . N with — oo < < oo and fixed values of r n , chosen such 
that all the singularities of the transform DID(Sq, h, q) are coordinate- wise on the left 
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of the lines C n . Many algorithms approximate the integrals in (|23p by a finite linear 
combination of the transform at some specific nodes with certain weights. Three 
approaches have been studied by Abate et al. [2], based on Fourier series expansion, 
combinations of Gaver functionals, and deformation of the integral contour. Here 
we concentrate on the last method developed by Talbot [31J, since reports in the 
literature (e.g. [2]) suggest that this approach offers high performance for a short 
time of execution, which our numerical results confirm. We write 

DID(S Q , h, T) = J* ■ ■ ■ J* p x {0) ■ ■ ■ (3 N (9)DID(S , h, q(0))dfl, (24) 

with n = 1, ... TV, (3 n (9) = w n e irnWnTn , q n {9) = ir n w n , and 

w n = -1 + iO + i(0 cot 6-1) cot 9. 

Since DID is a real valued function, DID is also equal to the real part of the integral 
on the right-hand side of (|24|) . which can be used to reduce the calculation by a 
factor of two. To illustrate the evaluation of the integrals (|24p . we present concrete 
expressions for the approximating sums when iV = 4 (which is the setting that will 
be implemented later on). Defining 



we obtain 



9 k n = kn/M and 



M-l M-l M-l M-l 



2M 

52V 



DID(S ,h,T) « ^f^fr E E E E 

J * fci=0 k 2 =0 k 3 =0 ki=0 

l3k 1 /3k 2 f3k 3 (3k 4 f{qk 1 /T 1 ,q k2 /T 2 ,qk 3 /T 3 , q k jT 4 ) + ]3 kl f3k 2 (3k 3 l3k i f(q kl /T 1 ,q k2 /T 2 ,q k jT 3 , q ki /T 4 ) 
+ l3k 1 P k2 ^ k3 /3 k4 f{qk 1 /T 1 ,q k2 /T 2 ,q k jT 3 ,q ki /T 4 ) + li kl P k J3 k3 P k J{q kl /T 1 ,q k2 /T 2 ,q k jT 3n q kl /T 4 ) 
+ Pk 1 Pk2pk 3 PkJ{lh 1 /Ti,qk 2 /T 2 ,qh 3 /T 3 ,q ki /T 4 ) + ]3 k j3 k2 (3 k3 l3 k J(q ki /T 1 ,q k2 /T 2 ,q k jT 3 ,q k4 /T 4 ) 
+ PkJl k2 fik 3 ll k J(qkjT\Ak 2 l T ^qk 3 lT 3l q k jT 4 )+J k ]i k ^ 

where / is equal to DID. The weights and the nodes are given by 

2M 2kir 
qo = -jp Qk = — (cot(/c7r/M) + i), < k < M, 

(3 = 0.5e 90 , fa = (1 + i(fc7r/M)(l + [cot(A;7r/M)] 2 ) - i cot(/c7r/M))e (?fc . 

Since the weights and nodes are independent of the transform, the calculation time 
of the algorithm can be reduced by pre computing and storing weights and nodes. 
The speed of convergence and the accuracy of the Talbot algorithm will depend on 
the regularity of the Laplace transform /. Although universal error bounds are not 
known, Abate et al. pQ showed numerically that the single parameter M can be used 
to control the error and can be seen as a measure for the precision. They found 
after extensive numerical experiments that for a large class of Laplace transforms 
the relative error is approximately 10 -0 ' 6M . For high dimensional inversion, extra 
accuracy in the inner sums may be needed to obtain a sufficient degree of precision 
for the outer sums, which can be achieved by increasing M. 
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C.2 Fractional Fourier Transform 



To evaluate down-and-in call option prices (DIC), we invert the Fourier-Laplace trans- 
form ()18|) over log-strike and time periods. For the inversion of the Laplace transform 
we again apply the Talbot algorithm. In the case of two time periods, with 

/o(?i)?2) = DIC (S ,h,v,q), 

we find that the Fourier transform DIC* can be approximated by the following sums: 

DIC*{S ,h,v,T) 

j M-lM-l 

~ 5 2 TT Y {PkiPtoMqkjTu q k2 /T 2 ) +P kl P k2 f v (q k jTi,q k2 /T 2 ) 

1 2 fc 1= 0fc 2 =0 

+ Pk 1 h 2 fv{qkjT 1 ,q k jT 2 ) + /3 kl ]3 k J v (q kl /T 1 ,q k2 /T2)} . 

Unlike the case of the inversion of DID, we cannot reduce the calculation time by 
two by using complex conjugates, since the function DIC* is not real valued. Down- 
and-in call prices are then obtained by inverting the Fourier transform over strike: 

„— ak roo 

DIC (S , h, k, T) = / e- ivk DIC*(v)dv, 

n Jo 

where a is the rate of exponential dampening. This integral is approximated for a 
set of log-strikes between (— xo,xo) as a summation: 

q -ak -rT N ~ 1 

DIC(S , h, k, T) « — V w j e-W- X0+kX )DIC*(5j)6, k = 1 • • • N - 1, 

IT ' ^ 

(25) 



7T 

3=0 



where (wj)j~ are the integration weights defined by the trapezoidal rule with wq = 
wn-i = 0.5 and 1 otherwise, A = 2xq/N is the log-strike grid step-size and 5 is the 
u-grid step-size. Carr and Madan [9J and Chourdakis [13] set 5 = 0.25. 

To have accurate prices for any strike, the log-strike grid spacing A needs to be 
sufficiently small. A common approach is to apply directly the Fast Fourier Transform 
(FFT) and to compute the summation ([25]) on a fixed log-strike range (— xo,xo) 
with xq = tt/5 using many points N. Bailey and Swarztrauber [5], [B] propose 
an alternative approach, and define the Fractional Fast Fourier transform (FrFFT), 
which uses an arbitrary range. Chourdakis [13] showed that the FrFFT can be used 
to calculate option prices with less points without losing accuracy. He reported that 
the FrFFT is 45 times faster than the FFT for the calculation of European option 
prices. Since in our case the Fourier transform DIC* is obtained numerically, we 
chose to employ the FrFFT. We now briefly specify the form of this algorithm in our 
setting, and refer for further details to [5], [6], [13]. The resulting sum is then given 
by 

q -(ak+mk 2 u) -rT N ~ 1 

DIC(S ,h,k,T) « 1 y w je -^ !2 " e -*Kk-j) 2 "DIC*(5j)5, 

j=0 
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where k = 1 ■ ■ • N — 1, Wj = Wje lx ° s i and v = 5xq/Nit. Extending this summation 
into a circular convolution over 2N yields 



C e -(ak+ink 2 u) -rT 2JV-1 

DIC{S ,h,k,T) « -5 ^VjZk-j, k = l---N-l, 

71 3=0 



where 

= {vje-^ 2v DIC*(Sj)6, Zj = e-^ 2 ", j < N, 

and 

y,=0, z,=e-^- 2 ^, j>iV. 

This equation can be rewritten in terms of three discrete Fourier transforms: 

Q ne -{ak+mk 2 v) -rT 
DIC(k) « ^ _ ^ F -i (Ffc(2/)Ffe(2))) fc = l--iV-l, 

with 

JV-l N-l 

F{x) = Xje-W/", F~\x) = x^ k ' N . 

3=0 3=0 

Although the latter sum is computed by invoking two Fourier transforms and one 
inverse Fourier transform, this approach has the advantage of computing the option 
prices on a specific log-strike window (— xq : x$) with independent grids S and A and 
requires less points. 
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